Formation of density singularities in ideal hydrodynamics of freely cooling inelastic 

gases: a family of exact solutions 
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We employ granular hydrodynamics to investigate a paradigmatic problem of clustering of par- 
ticles in a freely cooling dilute granular gas. We consider large-scale hydrodynamic motions where 
the viscosity and heat conduction can be neglected, and one arrives at the equations of ideal gas 
dynamics with an additional term describing bulk energy losses due to inelastic collisions. We em- 
ploy Lagrangian coordinates and derive a broad family of exact non-stationary analytical solutions 
that depend only on one spatial coordinate. These solutions exhibit a new type of singularity, where 
the gas density blows up in a finite time when starting from smooth initial conditions. The density 
blowups signal formation of close-packed clusters of particles. As the density blow-up time t c is 
approached, the maximum density exhibits a power law ~ (t c — t)~ 2 . The velocity gradient blows 
up as ~ — (t c — t)^ 1 while the velocity itself remains continuous and develops a cusp (rather than 
a shock discontinuity) at the singularity. The gas temperature vanishes at the singularity, and the 
singularity follows the isobaric scenario: the gas pressure remains finite and approximately uniform 
in space and constant in time close to the singularity. An additional exact solution shows that the 
density blowup, of the same type, may coexist with an "ordinary" shock, at which the hydrody- 
namic fields are discontinuous but finite. We confirm stability of the exact solutions with respect 
to small one-dimensional perturbations by solving the ideal hydrodynamic equations numerically. 
Furthermore, numerical solutions show that the local features of the density blowup hold universally, 
independently of details of the initial and boundary conditions. 

PACS numbers: 45.70.Qj, 47.40.-x 



T3 

C 

O 

i o i 

> 
00 

o 

O 

o 
o 



X 



I. INTRODUCTION 

Structure formation in many-body systems is one of 
central problems of non-equilibrium physics. A most 
spectacular phenomena of structure formation is clus- 
tering of matter. Here an initially structureless, al- 
most homogeneous distribution of particles of matter 
self- assembles into clusters. Such an evolution cannot 
proceed indefinitely in a gas where interactions between 
the particles are (i) short-range and (ii) Hamiltonian. 
There are two important classes of gases where one of 
these two properties is violated, and clustering can oc- 
cur. In the first one long-range forces, such as gravity, 
are present. Clustering provides here a natural mech- 
anism of star formation 1 - and of the large-scale struc- 
ture of the Universe^,. The second class of systems are 
dissipative systems where interactions between the par- 
ticles, at the level of an effective description, are non- 
Hamiltonian. One well-known example is optically thin 
gases and plasmas that cool by their own radiation, and 
dense condensations develop 3,4 ' 5 . In this paper we con- 
sider a more basic non-Hamiltonian many-body system: 
a gas of inelastically colliding macroscopic particles, or 
granular gas. Here particles lose energy to their internal 
degrees of freedom. 

The granular gas is the low-density limit of granular 
fiows££ In its simplest version, the granular gas model 
deals with a dilute assembly of identical hard spheres 
(with diameter a and unit mass) who lose energy at in- 
stantaneous binary collisions in such a way that the nor- 
mal component of the relative velocity of particles is re- 



duced by a constant factor < r < 1 (the coefficient of 
normal restitution) upon each collision. Granular gases 
exhibit a plethora of structure forming instabilities, in- 
cluding the clustering instability of a freely cooling ho- 
mogeneous inelastic ga^a^d^^^e^,!^ This 
instability brings about the generation of a macroscopic 
flow and formation of dense clusters of particles. 

The natural language for a theoretical description of 
macroscopic flows of a granular gas is the Navier-Stokcs 
granular hydrodynamics 6 ' 7 . Although the criteria of its 
validity are quite restrictive, see below, granular hydro- 
dynamics is instrumental for theoretical investigations, 
and often prediction, of a host of collective effects in 
granular flows. Recently, applications of granular hy- 
drodynamics have been extended to non-stationary flows 
of granular gase o 17 ' 20 i 21 ' 22 . The non-stationary settings 
provide sharp tests to continuum models of granular flows 
and help evaluate their domains of validity. This is es- 
pecially true when the time-dependent solutions of the 
continuum equations develop finite-time singularities^, 
such as the recently predicted density blowup in freely 
cooling granular gases: at zero gravit y 17 i 22 , and at finite 
gravity 21 -. 

We will assume throughout this paper nearly elastic 
particle collisions, a very small gas density (that we de- 
note by p), and a very small Knudsen number: 

l-r<l, po- d <rl, and i /ree /L<l. (1) 

Here d > 1 is the dimension of space, I free is the 
mean free path of the particles, and L is the charac- 
teristic length scale of the hydrodynamic fields. Under 
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these assumptions the Navier-Stokes hydrodynamics pro- 
vides a quantitatively accurate leading-order theory for 
granular gases^* 7 -. Previous work o 9 ' 10 employed hydro- 
dynamic equations to show that, for sufficiently large 
systems, the homogeneously cooling state of the gran- 
ular gas is linearly unstable. The unstable perturbations 
grow via two (linearly) independent modes: the shear 
mode corresponding to the development of a macroscopic 
solenoidal flow, and the clustering mode corresponding 
to the development of a potential flow that causes for- 
mation of clusters of particles. How does the cluster- 
ing mode develop beyond the linear regime, and do the 
hydrodynamic nonlincarities arrest the density growth? 
For the clustering instability in large systems (as well as 
for the gravitational instability 2 ) the growing perturba- 
tions bring the system into a fully developed non-linear 
regime & n ' 14 i 16 ' 17 i 18 . Therefore, one has to deal with fully 
non-linear hydrodynamic equations. Not restricting our- 
selves to a close proximity to the homogeneously cool- 
ing state, we can formulate the problem in a more gen- 
eral way and explore different nonlinear flows of a freely 
evolving granular gas. 

Solving the hydrodynamic equations analytically, and 
even numerically, is a difficult task, and additional simpli- 
fications are needed. We will make two additional simpli- 
fying assumptions in this work. First, assuming a large- 
scale flow, we will neglect the viscous and heat conduc- 
tion terms in the hydrodynamic equations. Second, we 
will assume that the macroscopic flow is one-dimensional. 
A natural environment for the second assumption is pro- 
vided by the geometry of a narrow channel with perfectly 
elastic side walls that we adopt here, following Refs . 17 ' 18 . 
Although the microscopic motion of particles in the chan- 
nel remains two- or three-dimensional (2d or 3d), the 
macroscopic flow depends only on one spatial coordinate: 
the coordinate along the channel. This is because, in a 
narrow channel, both the shear mode and the cluster- 
ing mode in the transverse directions are suppressed (see 
Refs. [T7IIT8I for detail). As a result, one can focus on the 
development of the (one-dimensional) clustering mode as 
it enters a strongly nonlinear regime. 

Working in the channel geometry, Efrati et alM con- 
sidered the long-wavelength limit of the clustering insta- 
bility, when the linear growth rate of the instability is 
the highest. In this case the inelastic energy loss of the 
gas is the fastest process, and the gas pressure drops, al- 
most instantaneously, to a very small value. The further 
dynamics are then (almost) purely inertial which would 
bring about a finite-time blow-up of the velocity gradi- 
ent and, therefore, of the density. This is the well known 
blow-up of the free flow 2 ^. Its signatures were observed in 
the numerical solution of the hydrodynamic equations 17 
until the maximum gas density became so high that the 
numerical scheme lost accuracy. The numerical results of 
Ref.— were tested in molecular dynamics (MD) simula- 
tions in 2d of a freely cooling gas of inelastically collid- 
ing disks in a long and narrow channel. The MD sim- 
ulations supported the free-flow blow-up scenario until 



the time when the gas density approached the hexagonal 
close-packing value, and the further density growth was 
arrested. 

As we have recently found^ 2 -, the free flow asymptotics 
does not hold all the way to the density blowup. Very 
close, in time and in space, to the (attempted) free-flow 
singularity, the compressional heating starts to act and 
temporarily stabilizes the gas temperature at a small 
but finite value. As a result, the gas pressure again be- 
comes important: it breaks the purely inertial dynamics 
and, though unable to stop the density blowup, dramat- 
ically changes the local blowup properties. It turns out 
that, after a brief crossover, the further dynamics obey a 
new blowup scenario, not limited to the long-wavelength 
limit 2 ^. The new scenario has the following features. As 
the blow-up time t c is approached, the maximum den- 
sity exhibits a power law ~ (t c — t)~ 2 . The velocity 
gradient blows up as ~ — (t c — t) _1 , whereas the veloc- 
ity itself remains continuous and develops a cusp, rather 
than a shock discontinuity, at the singularity. The gas 
temperature vanishes at the singularity, but the pressure 
there remains finite. This blowup, which obeys the above 
asymptotic laws near the singularity, emerges universally, 
that is for generic initial and boundary conditions. (Note 
that, in the long-wavelength limit, the crossover from the 
free flow regime to the finite-pressure regime does not 
happen when the initial gas density is not very small. 
In this case excluded particle volume effects interfere in 
the dynamics, and arrest the density growth, before the 
crossover has a chance to occur, as indeed was observed 
in Ref^.) 

These findings are based on extensive numerical sim- 
ulations (numerical solutions of the hydrodynamic equa- 
tions for a host of initial and boundary conditions) and a 
family of exact analytical solutions of the hydrodynamic 
equations, without and with shocks, briefly announced in 
Ref^ 2 -. In the present work we give a detailed account of 
the exact solutions and investigate their structure close 
to the singularity. We verify the exact solutions numer- 
ically. We show that the singularity follows the isobaric 
scenario: the gas pressure is approximately uniform in 
space and constant in time in a close vicinity of the de- 
veloping singularity. Finally, we evaluate the validity of 
the exact solutions close to the singularity by estimating 
the role of additional physical processes: the viscosity, 
heat conduction and excluded particle volume effects. 

The remainder of the paper is organized as follows. In 
Section II we introduce the equations of ideal granular 
hydrodynamics (IGHD) and discuss their general prop- 
erties. In Section IlIII we employ Lagrangian coordinates 
and derive a family of exact solutions of the IGHD equa- 
tions, consider some particular cases of the solutions and 
investigate global and local properties of the solutions. 
In Section [TV] we adapt the exact solutions to describe a 
different setting, where a piston moves into a granular gas 
at rest, and a density blowup, developing at the piston, 
coexists with an "ordinary" shock wave propagating into 
the gas. Besides demonstrating the presence of the two 
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different types of singularities in the same system, this so- 
lution allows an arbitrary initial density profile at large 
distances, showing that the density blowup is a local pro- 
cess. Section|V]presents the results of numerical solutions 
of the IGHD equations that confirm stability of the exact 
solutions with respect to small one-dimensional pertur- 
bations and establish universality of the density blowup 
for different initial conditions. In Section [VII we discuss 
the role of non-ideal effects, neglected in our solutions, 
close to the singularity. In Section [VlT] we summarize our 
results and discuss their bearing on cluster formation. 



II. IDEAL HYDRODYNAMICS OF A FREELY 
COOLING GRANULAR GAS 

Under the three strong inequalities {T]), the Navier- 
Stokes granular hydrodynamics provides a quantitatively 
accurate leading-order theory for granular gases. It deals 
with three coarse-grained fields: the mass density p(x, t), 
the mean flow velocity v(x, t) and the granular tempera- 
ture T(x, t). An additional coarse-grained field, the pres- 
sure p(x, t) , is related to the density and temperature 
by the perfect gas equation of state p = pT. Assuming 
a one-dimensional macroscopic flow, we can write these 
equations as 



dp d{pv) _ 
dt dx 
' dv dv 



0. 



p ^t +v Tx 



d{ P T) 
dx 



(2) 



dT 
~dt 



dT 

dx 



kq d 
p dx 



-( 7 -l)T 



dT 

dx 



dv 



A P r 3 / 2 

dx 

1/0(7- l)Vr fdv 
p \dx 



(4) 



Here 7 is the adiabatic index of the gas (7 = 2 and 5/3 
for d = 2 and d = 3, respectively), A = 2it {d - 1 V 2 (l - 
r 2 )<T d ~ 1 /[dT(d/2)} (see e.gM), and T(. . .) is the gamma 
function. In addition, v$ = {2a ^/n) -1 and kq = 4vq in 
2d, and vq — 5(3cr 2 Y / 7r) _1 and kq — 15fo/8 in 3d, see 
Ref. @. The only difference between Eqs. ©-(H]) and 
the standard gas dynamic equations for a dilute gas of 
elastically colliding spheres is the presence in Eq. ([4]) of 
the inelastic energy loss term — ApT 3 / 2 . 

There are three types of dissipative terms in Eqs. ([2])- 
(j4|): the viscous terms in Eqs. |3|) and (J4j) , the heat 
conduction term in Eq. (0} and the energy loss term in 
Eq. ((3]). The viscous and heat conduction terms include 
spatial gradients of the hydrodynamic fields, whereas 
the energy loss term is independent of the gradients. 
When the characteristic hydrodynamic length scale of the 
flow is sufficiently large, the viscous and heat conduction 
terms can be neglected, while the energy loss term should 
be kept, and we arrive at the equations of ideal granular 
hydrodynamics (IGHD): 

dp . d(pv) fdv dv\ d{pT) 

di + ^x- = °> P {-di +V d-x )=-^' (5) 



dT dT , , dv k 

m +v ^ = -^-^ T 8- x - ApT/ - 



(6) 



For consistency, all the assumptions must be checked a 
posteriori, after a hydrodynamic problem in question is 
solved, and the hydrodynamic length and time scales are 
found. 

The basic state of the freely cooling gas is the homo- 
geneously cooling state described by the Haff law 25 : 



P = Pa, v = 0, T = 



To 



1 + A Po T Q 1/2 t/2 



(7) 



This state corresponds to the initial conditions p(x, t — 
0) = po = const, T(x, t = 0) = Tq = const and v(x, t) — 
0. Obviously, the IGHD equations reproduce the Haff's 
law exactly, as the homogeneously cooling state does not 
include gradients of the hydrodynamic fields. 

A more meaningful example of a situation where the 
IGHD applies is provided by the linear stability anal- 
ysis of the homogeneously cooling state. This analy- 
sis, in the framework of the complete, non-ideal Navier- 
Stokes hydrodynamics ©-© and its extensions was per- 
formed by many workers, starting from Goldhirsch and 
collaborators^ and McNamara— . The main results of the 
linear stability analysis can be described, at the level of 
order-of-magnitude estimates, as follows. The evolution 
of a small sinusoidal perturbation with the wave num- 
ber k is determined by two competing processes. The 
inelastic energy loss tends to enhance the fluctuations 
on the cooling time scale [(1 — r 2 )cr d_1 pVT}^ 1 which 
is fc-independent. In its turn, the viscosity and ther- 
mal conduction tend to erase the perturbation on a time 
scale [k 2 l 2 free po- d - 1 VT]- 1 , where l free ~ l/ina^ 1 ) is 
the mean free path and l 2 ' ree po~ d ~ 1 VT is the character- 
istic value of the viscosity/hcat conductivity. Balanc- 
ing these two time scales, one obtains the critical wave 
number k c ~ Ijree^^ ~ r so that the perturbations with 
k < k c grow, while the short- wavelength perturbations, 
k > k c , decay^iiS. In a narrow channel, such that per- 
turbations in the transverse direction with wave num- 
bers smaller than k c are not allowed, only perturbations 
along the channel will grow. As a result, by the end of 
the linear stage, the hydrodynamic fields are effectively 
one-dimensional and have characteristic length scales of 
the order of lf ree /y/l — r or longer. The validity of the 
Navier-Stokes hydrodynamics ©-((31) f° r the description 
of the whole range of unstable perturbations demands 
a strong inequality \/l — r <C 1. On the other hand, 
the IGHD model ([5]) and © is valid if we demand a 
strong inequality k <^ k c . Indeed, one can check that, in 
this case, the growth rate of the clustering instability, as 
obtained from the IGHD, coincides in the leading order 
in k/k c with that obtained from the full hydrodynamic 
equations. 

In Section |VT1 we will perform a consistency check, and 
establish the validity domain, of our exact nonlinear so- 
lutions of the IGHD equations ([5]) and ©. Now let us 
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discuss the basic properties of these equations. Although 
much simpler than the non-ideal equations ©-(HI), the 
nonlinear IGHD equations still present a hard mathe- 
matical problem. Going back to elastic particle collisions, 
A = 0, one recovers the ordinary ideal gas dynamics in 
one dimension. Among most interesting solutions are 
those describing the development of wave-breaking singu- 
larities when starting from smooth initial data 2 ^ 2 - 7 . Note 
that, even for A = 0, the general initial value problem is 
not soluble analytically, except for the particular case of 
an isentropic flow, where the entropy per unit volume 
s(p, T) — pin (T//9 7-1 ) is uniform in space and constant 
in time2&2Z. Needless to say, at A > 0, Eqs. © and © 
do not allow isentropic solutions. The total entropy of 
the fluid S = J s(p,T) dx, governed by Eq. © and ©, 
is monotone decreasing: 



dS 
dt 



^ = -A / p 2 T 1 / 2 dx 



(8) 



where we have assumed that there is no net entropy 
flux from the boundaries. As expected from the micro- 
scopic picture, the local entropy loss rate in Eq. © is 
proportional to the particle collision rate. The entropy 
loss implies that the system may exhibit self-organization 
phenomena^ as is indeed observed in the process of clus- 
tering instability. 

Before deriving a family of exact solutions exhibit- 
ing a finite-time density blowup, we note two rescaling 
symmetries of Eqs. © and ©. The first symmetry 
relates solutions at different A > 0. If p(x,t), v(x,t) 
and T(x,t) solve Eqs. © and © at some A, then the 
rescaled fields p[(A'/A)x,(A'/A)t], v[(A'/A)x, (A'/A)i] 
and T[(A'/A)a;, (A'/A)t] solve the same system with the 
cooling coefficient A'. Therefore, the general mathemat- 
ical properties of Eqs. © and ©, such as the existence 
of singularities, are the same for any A > 0. The sec- 
ond symmetry relates solutions of Eqs. © and © at the 
same A. If p, v and T solve Eqs. © and ©, then p, v 
and T defined by 

p(x, t) = ap(ax, ot\f~f3t), v(x, t) = y/flv(ax, a.\/ ~/3t), 

f(x,t) =PT(ax, ay/fit), (9) 

with any a > and {3 > 0, also solve Eqs. © and ©. 
These symmetries are exploited in the following. 



III. DEVELOPMENT OF DENSITY 
SINGULARITIES 

A. Lagrangian coordinates and exact solutions 

Let us rewrite the governing equations © and © in 
terms of the pressure p — pT, rather than temperature: 



dp 9(pv) 



dt 
dp 
dt 



dx 
dp 
dx 



= 0, 



-IP 



' dv dv\ dp . n . 



dt 



dv 

O.r 



dx 

1/2 3/2 



dx ' 



Ap'i'p 



(11) 



The family of exact solutions, presented in this Section, 
are smooth initially but become singular in a finite time 
t c . At the singularity, the density blows up, in contrast 
to the ordinary "wave-breaking" singularity of the ideal 
gas dynamics (A = 0), where only the gradients of the 
hydrodynamic fields blow up 26 ' 27 . An exact solution of a 
different type, presented in Section \W\ includes a shock 
already at t = 0. That solution also exhibits a density 
blowup, and the local properties of the blowup are the 
same as in the initially smooth solutions. 

Let us introduce Lagrangian coordinates. The coordi- 
nates x(m, t) of the fluid particles obey the equation 



dx(m, t) 
dt 



v(x(m, t),t), 



(12) 



where m is a continuous label (a Lagrangian coordinate) 
of particles. The defining property of the exact solutions 
that we are going to present is independence of the par- 
ticle accelerations dfx(m,t) of time: the fluid particle 
coordinates x(m, t) satisfy the equation 



x(m, t) = x(m, 0) + v(m, 0)t 



a{m,0)t 2 



(13) 



where v(m, 0) and a(m,0) are the initial velocities and 
accelerations of the fluid particles, respectively. As we 
will see later, the total mass of gas is finite for our solu- 
tions. As a result, the pressure must vanish at the bound- 
aries, implying existence of a point with a zero pressure 
gradient (and, therefore, a zero particle acceleration) in 
between. A fluid particle that has a zero acceleration 
(that is conserved in our solutions) moves with a con- 
stant velocity, and we choose to work in such a frame of 
reference where this particle is at rest, so that the condi- 
tions v(x = 0, t) = Q and d x p(x = 0, t) = are obeyed. 

It is convenient to choose the Lagrangian coordinate 
to to be the mass coordinate 2 ^: 



m(x, t) 



p{x' , t)dx 



(14) 



that is the mass content between the Eulerian points 
and x. The inverse transformation x(m, t) is 



x{m, t) 



m dm' 



o p{m',t) 



(15) 



In the Lagrangian coordinates Eqs. (fTTJ|) and (TiTj) look 
simpler: 



d ( 1 \ dv dv dp 

dt \p J dm ' dt dm ' 

9v 3/2 1/2 
dt am 

Let us calculate the Lagrangian acceleration: 
d 2 x{m,t) dp(m,t) 



dt 2 



dm 



(16) 
(17) 

(18) 
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where we have used Eqs. (|12j) and (|16|) . For the solu- 
tions obeying Eq. (fT3"|) , the Lagrangian pressure gradient 
d m p(m, t) should be independent of time. Since the pres- 
sure p(m, t) is time-independent (and zero) at the gas 
boundaries, it can depend only on to. Then Eqs. (fT6|) 
and |n]) yield 



9 2 p 
dm 2 



= —p 2 p, where p = 



A 



7 



V2- 



(19) 



In view of the zero acceleration at the origin, d m p(0, t) = 
0, Eq. (fT9|) yields p = 2Acos(pm), where A is constant. 
The resulting family of solutions is 



p(m,t) 
p(m,t) 

v(m, t) 



2 A cos(pm), 

p(m,0) 



[1 — pt^J Ap(m, 0) cos(^to)] 5 
2A^ii sin(/im) , 



(20) 
(21) 



(22) 



where we have used the condition v(m = 0,t) = 0. This 
family of solutions describes a compression flow (the ve- 
locity gradient is negative everywhere), as the compres- 
sional heating is balanced, in Eq. (fT7|) . by the inelastic 
cooling. The solutions include an arbitrary non-negative 
function p(m, 0): the initial gas density. The arbitrary 
constant A > that, together with p(m,0), determines 
the (time-dependent) Mach number of the flow, appears 
due to the rescaling symmetry of the equations and corre- 
sponds to the constant (3 in Eq. ([9]). One can check that 
the constant a in Eq. © corresponds to the freedom of 
multiplying p(m, 0) and A by a. 

As the pressure p must be non-negative, and vanish 
at the (finite or infinite) boundaries of the freely moving 
gas, the solutions (|2H)l can hold only on a finite interval 
(— n/2p, Tr/2p) (we assume that the gas region is single- 
connected, and the interval includes m = 0). Therefore, 
the total mass of the gas in these solutions is finite and 
fixed by parameters A and 7: p{x, 0)dx — n/p — 

V2n-//A. 

Once the solutions (f2"0"l) - ([2"2"]) in the Lagrangian coor- 
dinates are known, we can return to the Eulerian coor- 
dinates by using, at any time t, Eq. (fT5|) . Depending 
on the particular choice of the initial density, there are 
two possible types of solutions (f2"0")) - (|2"2"|) . First, the fixed 
mass of the gas ir/p can be distributed, at t = 0, over 
either an infinite, or a finite x-interval. This is deter- 
mined by the behavior of p(m, 0) near m = ± ir/2p. For 
example, let p(m,Q) ~ (tt/2/x — m) 1+a near to = w/2p. 



Then it follows from Eq. (|T5|) that the gas occupies an 
infinite (correspondingly, a finite) interval of positive x 
if a > (correspondingly, a < 0). The velocity can 
be either finite, or infinite at the gas boundaries. For 
example, for the same behavior of the initial density 
p(m,t = 0) ~ (t/2// — m) 1+a near m = n/2p one ob- 
tains a finite (correspondingly, infinite) gas velocity at 
to = 7f/2/i for a < 2 (correspondingly, a > 2). 



Now let us consider what types of initial conditions 
evolve according to Eqs. (|2T)1) - ({2"2"1) and discuss the density 
blowup that is brought by this evolution. 



B. Initial conditions and density blowup for the 
exact solutions 



A particular member of our family of exact solutions 
(f20)) - (f2"2"l) is determined by a specific choice of the con- 
stant A > and of the initial density p(m, 0) > 0, de- 
fined on the interval [— ir/(2p), ir/(2p)]. In the Eulerian 
coordinates one can specify an arbitrary initial density 
profile p(x, 0) that has a fixed total mass n/p: 



p(x, 0)dx = n/p. 



(23) 



Once p{x, 0) and A are specified, one needs to choose the 
origin so that the gas masses to the left and to the right 
of the origin are the same [and equal to 7t/(2p)]. Then 
the initial gas pressure in the Eulerian coordinates is 



p{x,0) =2Acos[ p p(x',0)dx' 



(24) 



Now, using Eq. (|2"2"]) , we calculate the velocity gradient 
in the Eulerian coordinate at t = 0: 



dv(x,0) dv(m,0) 
= P{m,0)- 



dx 



dm 



= —2py / Ap(m, 0) cos(pm). 



(25) 



which, in view of the condition v(x = 0, t) — 0, yields 



v{x,0) = -p v / 2p(x',0)p{x',0)dx' (26) 
^0 



withp(x,0) from Eq. Ipi|) . Equations flM} and (J26l) show 
that, once p{x, 0) is a smooth function of x, then the 
initial pressure and velocity are smooth functions as well. 

Let us now proceed to the properties of the solutions. 
As we already noted, these solutions describe a motion of 
fluid particles with a time-independent acceleration, see 
Eq. (fl3)) . This time- independent acceleration is 



a(m, 0) 



9p(m, 0) 
dm 



= 2pA sin(/iTO,). 



(27) 



The evolution described by Eqs. ([20[) - (f22]) brings about a 
singularity of this initially smooth flow: p(m, t) blows up 
in a finite time, while the rest of the flow variables - the 
pressure and velocity - remain finite. The density blow up 
occurs at the Lagrangian point too where p(m, 0) cos(pm) 
reaches its maximum: 

p(too, 0) cos(/iTOo) = max [p(m, 0) cos(pm)] . (28) 

Interestingly, the point too corresponds, in view of 
Eq. (f2l))) . to the point of the absolute (negative) mini- 
mum of the velocity gradient in the Eulerian coordinates, 
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just as in the case of the free flow (that is, zero pressure) 
singularity 24 . The singularity occurs when the Jacobian 
of the Lagrangian transformation of the fluid particles 
vanishes for the first time: d m x(m — mo,t c ) = 0. For 
the time t c and the Eulerian coordinate of the singularity 
xq we find 



tc 



1 



p^/Ap(m a ,0) cos(/im ) ' 
. f mo dm 1 
Jo p{m',t c ) 



(29) 



At the (fixed) Lagrangian point of singularity mo the 
density blows up as 



p(m ,t) 



p(m o ,0) 
(1-i/tc) 2 



(30) 



In the Eulerian coordinates the blowup develops, in gen- 
eral, in a moving point: 



p[x(m ,t),t] = 



p[x(m , 0),0] 

d-t/t c y ■ 



(31) 



The velocity gradient at the singularity point x(mo,t) 
diverges. The divergence law can be easily found in the 
Langrangian coordinates, by using the continuity equa- 
tion and Eqs. flUI) and ([29]) : 



dv(x, t) 



dx 

d_ 

ot 



x—x(m,Q,t) 

[In p(m,t)} 



p(m,t) 



dv{m, t) 



dm 



tr.-t 



(32) 



Finally the (finite) pressure is conserved on the La- 
grangian trajectory: 



p(x = x(iriQ,t),t) = 2Acos(pmo). 



(33) 



Note that the pressure does not have a minimum at the 
singularity point, so this flow does not conform to the 
popular "pressure instability" scenario^. 

It is instructive to consider several particular examples 
of solutions starting with the mass distributed over an 
infinite interval of x. 



C. Solutions with mass distributed over an infinite 
interval of x 

In our first example the initial density profile in the 
Lagrangian coordinates is p(m, 0) = pocos(pm). To re- 
turn to the Eulerian coordinates, we use Eq. (|15p and 
obtain 



sinh 



x(m, 0) 



= tan(/im) 



(34) 



where we have introduced the characteristic inelastic 
cooling length scale I = l/(ppo) whose meaning will be- 
come clear shortly. We use Eq. (jM)) to express pirn, 0) = 



po cos(pm) through x. The rest of initial conditions fol- 
low from Eq. (J2UJ) and Eq. ^ at t = 0. Note that 
the initial gas temperature T(m, 0) = p(m, 0)/p(ra, 0) = 
2A/pq = To is uniform in space. The initial conditions 
are 



= coshCr/0' 
v(x,Q) = — \J 2Xo arcsin 



T(x,0)=T , 
t a nh(|) 



(35) 
(36) 



The initial velocity profile describes an inflow of gas 
from plus and minus infinity with a finite velocity there: 
lim-c^ioo v(x, 0) = =F TTy^To/2. Now let us introduce the 
characteristic inelastic cooling time 



lV2 V2 



27 



To ppoV% Po^y/To ' 



(37) 



see Eq. ([6]). In its turn, / is the characteristic length 
scale the particles pass during the time r while moving 
with thermal velocity. According to Eqs. (|20p - (f2"2")l , the 
hydrodynamic fields in the Lagrangian coordinates evolve 
in time in this example as 



p(m,t) 
p(m,t) 

i>(to, t) 



p T cos(p,m), 
po cos(^m) 



(38) 

[l-(</r)cos( M m)] 2 ' (39) 
— \pXT~o~ [pin — (t/ t ) sin(/iTO,)] , (40) 



as depicted in Fig. [TJ 

Using Eqs. (fl~5|) and (j39|) we find the law of motion 
([T3"|) of Lagrangian particles, 



x(m, t) 



1 



3in(/im) \ 



1 — sm(pm) J 

f t\ 2 
pm + I — ) sin(/im) 



which, combined with Eqs. (13"9"1) - (|4T)]) . yields a parametric 
representation of the solution in the Eulerian coordinates, 
see Fig.[5J The density singularity occurs at x = at time 
t = t: 



p(o,t) 



Po 



(1-t/r) 2 
p(0, t) = poTo — const, 



T(0,t) 
dv 



dx 



(0,i) 



To(l-</r) 2 , 
2 



T-r 



(41) 



while v(x = 0,t) = 0. Notice that the pressure has a 
local maximum at the density blowup point x = m = 0. 

The above solution can be immediately generalized. 
Indeed, it is a particular case of a one-parameter family 
of solutions generated by the initial density profile 



Po 
2 



cosh 



(^ + a)+cosh- 1 (^- a ) 



(42) 

where a > is an arbitrary parameter. For a < a cr = 
arcsinh(l) = 0.88137 . . . there is a single density peak at 
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FIG. 1: (Color online) An example of a single-peak solution 
in the Lagrangian coordinates [Eqs. (|38|l - (|40| l]. and a numer- 
ical solution for the same parameters. The analytical and 
numerical solutions are depicted by solid lines and symbols, 
respectively. Shown are the inverse density (a), velocity (b), 
pressure (c) and temperature (d) at times t = 0.62 (circles), 
0.77 (squares) and 0.99 (triangles) as functions of the La- 
grangian coordinate m. Shown in (e) is the inverse square 
root of the density at the singularity point m = (circles). 
Shown in (f) is the total energy of the gas versus time as 
described by Eq. (I60|l . The density, velocity, pressure and 
temperature are rescaled to pa, yfTo/2, poTo/2 and To/2, re- 
spectively. The Lagrangian mass coordinate is measured in 
units of M/n — 1/p. Time is measured in units of r, so the 
density blows up at t = 1. Details of the numerical solution 
are given in Section IVl 
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-0.1 X 0.1 -0.2 Ox 0.2 
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FIG. 2: (Color online) The exact solution from Fig. [T]in the 
Eulerian coordinates. Shown are: the density in the logarith- 
mic scale (a), the velocity (b), the pressure (c) and the tem- 
perature (d) at rescaled times t = 0.62 (the dashed-dotted 
line), 0.77 (the dashed line), 0.99 (the solid line) and the 
blowup time 1 (the dotted line in c) as functions of x. The 
x-coordinate is measured in units of I. The rest of units are 
the same as in Fig. [T] 



x = 0, while at a > a cr there are two symmetric density 
peaks at x — ± I arccosh(sinh a) . The density profile (jl2"|) 
obeys Eq. <[23|) : the total mass of the gas remains equal 
to 7r//i. Equations (fl~5|) and (jl2"|) yield 



sinh 



x(m, 0) 
I 



= coshatan(/xm). 



(43) 



By setting p(x, 0) = pqTq cos[/im(x, 0)] we obtain the ini- 
tial temperature 



T(x,0) 




sinh a 
cosh 2 (x/l) 



(44) 



while v(x, 0) can be found from Eq. ([2BJ) , Now we calcu- 
late the initial conditions in the Lagrangian coordinates. 
Using Eqs. and (O we obtain 



p(m, 0) = po cos(/im)^/ 1 — tanh 2 a cos 2 (pm). (45) 
Then Eq. dHJ) yields 

dm' 



v(rn, 0) 



(1 — tanh 2 a cos 2 to') 1 / 4 



(46) 



Though this integral can be expressed via the Appell hy- 
pergeometric function of two variables, the integral form 
is more visual. The time history of the hydrodynamic 
fields in the Lagrangian coordinates is shown in Fig. [3] 
To go over to the Eulerian coordinates, we calculate the 
law of motion (jT3J) of the Lagrangian particles: 



x(m, t) 



In 



2f | 



coshatan(/iTO,) + \J 1 + cosh 2 <ztan 2 (/iTO.) 

» m dm/_ 

ii (1 — tanh 2 a cos 2 m') 1 / 4 



— ) sin(pm), 



and use it together with the Lagrangian solutions, see 
Fig. m 

The development of singularity in this case depends 
on the parameter a. For a < a cr the density and the 
pressure peaks remain at x — at all times until the 
singularity, while the singularity is of the same type as 
that observed for a = 0: 

A) 



p(Q,t) = — 

(V cosh a — t/r) 2 

T(0, t) = T (V cosh a - t/r) 2 , 

dv _ 2 
dx 



P(0,t) = PqT , 



rvcosha — t 



(47) 



The density blowup time is t c — r\/cosha. Now let 
us consider the case of a > a cr with two symmet- 
ric off-center density peaks at t = 0. Interestingly, at 
a cr < a < a cr , where a cr = arcsinh(\/2) = 1.14621 . . ., 
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10' 



m 



FIG. 3: (Color online) An example of a two-peak solution 
in the Lagrangian coordinates. Shown are the exact solu- 
tions (f20)) - ((22|) with the initial density from Eq. (|45)) : the 
inverse density (a), velocity (b), pressure (c) and tempera- 
ture (d) at the rescaled times t = 1.1 (the dashed-dotted 
line), 1.3 (the dashed line) and 1.45 (the solid line) versus the 
Lagrangian coordinate m. The density blows up at the La- 
grangian points iarccos ^^/2/3 coth 1.5^ = ±0.44628 ... at 

the rescaled time t c = 3 3/4 2" 1/2 tanh 1.5 = 1.45896 . . .. The 
units are the same as in Fig. [1] 



the singularity still develops at x = 0. Indeed, it is the 
maximum of the function p(m, 0) cos(/xm), rather than 
of p(m, 0) that determines, in view of Eq. ([28]). the sin- 
gularity point. For a cr < a < a cr this maximum is still 
at m = 771q = Q. As time progresses, the two symmetric 
density peaks move toward the origin, reaching it pre- 
cisely at the time of singularity. The pressure still has a 
maximum at x = 0, while in view of Eqs. (|47|) the time of 



singularity is still t c — rvcosha. Instead of looking for 
the maxima of p(m), it is convenient to look for the min- 
ima of the inverse density which, by virtue of Eq. (|21[) . 
can be written as 



Po 



1 



p{ m it) cos(fim)\J 1 — tanh 2 a cos 2 (pm) 

x- -i +f£\ cos(Mm). (48) 

[1 — tanh a cos 2 (pm) J V / 

Differentiating this function with respect to m one can 
verify that at t — t\J cosh a the minimum is indeed at 
m = 0. Finally, at a > a cr the density blows up sym- 
metrically at two Eulerian points x — ± x(m Q ,t c ) ^ 0, 
where 

mo = (l/p) arccos ( y/2/3 cotha^j ^ 0. 



Using the first of Eq. ([29]) we find i c /r = 
33/4 2- 1 / 2 tanh a. In this case the pressure gradient is 
non-zero in the singularity point, so the pressure is nei- 
ther maximum, nor minimum. 
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FIG. 4: (Color online) The exact solution from Fig. [3] in the 
Eulerian coordinates. Shown are the density in the logarith- 
mic scale (a), velocity (b), pressure (c) and temperature (d) 
at rescaled times t = 1.1 (the dashed-dotted line), 1.3 (the 
dashed line) and 1.45 (the solid line) as functions of x. The 
insets in (a) and (d) zoom on the density and temperature 
profiles, respectively, at t — 1.45 and clearly show the pres- 
ence of two symmetric peaks at 1 / developing density 
blowups simultaneously. The units are the same as in Figs. [1] 
andd 



D. Solutions with mass distributed over a finite 
interval of x 

If p(m, 0) does not vanish at \m\ — 7r/(2/z), or vanishes 
slower than linearly there, then the integral in Eq. (fT5]) 
converges, and the total gas mass n/ p is distributed over 
a finite interval of x. It follows from Eq. (|22|) that the 
velocity (which vanishes at x — m = and has a negative 
gradient everywhere) takes finite values at the ends of the 
interval over which the mass is distributed. As a result, 
the x-interval, occupied by the gas, shrinks in the course 
of evolution. Assuming for simplicity that p(m, 0) is an 
even function of m, we find from Eqs. (fT5|) and (f2l|) that 
the interval [— L(t) , L(t)] , occupied by the gas, shrinks 
with time as 



L(t) = L(0)~2pt 



\ I n ' dm' + pAt 2 49 



and reaches its minimum at t = t c . This minimum is al- 
ways positive except in the degenerate case of p(m, 0) = 
Po/ cos(pm), when the whole gas collapses into the point 
x = at the time of singularity. It turns out that, in this 
degenerate case, the solution is self-similar in the Eule- 
rian coordinates x and t and separable in the Lagrangian 
coordinates m and t. It can be shown that all self-similar 
solutions with a finite energy have an infinite density at 
some locations already at t — 0. Such initial conditions 
do not correspond to a dilute gas, so they will not be 
considered here. 
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A generic example of the solution on a finite Eule- 
rian interval is provided by the uniform initial density 
p{m, 0) — po- In the Eulerian coordinates this choice 
corresponds to the constant initial density, p(x,0) = po, 
at the interval [-wl/2,nl/2]. The solutions (pfj] ) -(f2"2 ]) be- 
come 



p(m,t) 
p(m,t) 

v(m, t) 



p T cos(pm), 
Po 



[1- (t/r)Vcos(^m)] 2 



2E (? 



• sm(pm) 



(50) 
(51) 

,(52) 



where E (. . . | . . .) is the elliptic integral of the second 
kind. The relation for x(m,t) is the following: 



x(m, t) ft 
, = pm - 4 - 
I \ T 



E (^ 



— I sin(/im). 



Here the singularity occurs at x — at time t c = r 
whereas L(t c ) > 0. The local structure of singularity is 
the same as in the case of an infinite interval, see subsec- 
tion [rrrFl 



E. Energy decay for the exact solutions 

A useful global characteristics of the clustering process 
is provided by the evolution of the total kinetic energy of 
all particles versus time, E(t) > 0. This quantity is con- 
venient to follow in experiment and in MD simulations. 
In the language of hydrodynamics it is 



E(t) 



7-1 



pv 
~2 



dx . 



(53) 



where the first term under the integral is the thermal 
energy density, and the second term is the macroscopic 
kinetic energy density. Let us compute E(t) for the exact 
solutions Eqs. (|20 |) - (f22)) . First, consider the conditions 
under which the initial total energy E(0) is finite. For 
the initial thermal energy we have 



E th (0) = f 

J —i 



= 2A 



p(x, 0)T(x, 0)dx 
"■/^ cos{pm)dm 



(54) 



tt/2 P (7- l)p(m,0) 



Whether this quantity is finite or not depends on the 
behavior of p(m, 0) near m = ±n/2p. For example, as- 
suming as before that p(m, 0) ~ {^/2p — m) 1+a near 
m = 7r/2/i, we find that E t h(0) is finite at a < 1 and 
infinite otherwise. A simple example of the initial con- 
dition with an infinite energy is p(m,0) = pocos 2 (pm) 
corresponding to the Lorentzian initial density profile 
p(x,0) = po/[l+(x/l) 2 } in the Eulerian coordinates. Here 



p(x,0)T(x,0) = poTo/Vl + (x/l) 2 



-00 < x < 00 . 



which is not integrable. Now, since at a < 2 the gas 
velocity is finite at m = ir/2p (see subsection llll A[) . then 
at a < 1 the initial macroscopic kinetic energy 



■Efcin(O) 



00 9 

pv 



dx = 



* /2f * v 2 (m,Q) 



TT/2fl 



dm (55) 



is also finite. Therefore, we assume a < 1 so that E(0) < 
00. Using the divergence form of the energy equation, 



d ( P T 



dt V7- 1 



pv 
~2 



d f jpvT 
dx \7 — 1 

A/0 2 T 3/2 



pv" 



7 -l ' 

we can express the decay rate of the total energy as 
A 



(56) 



dE 



7' 



1 



A 



p 2 T 3 / 2 dx 

-OO 

-ir/2p, P 



1/2 



(57) 



Using Eqs. $2U§ and |[2"T]). we obtain 



dE 
~dt 



A(2A) 



3/2 



7' 



1 



"P" dmcos 3 / 2 (pm) 
-TT/2 t i ^/p(m,0) 



At 



(58) 

It is easy to see that the integral in this equation con- 
verges for the assumed behavior of the initial density. 
Therefore, the energy decays quadratically in time: 



E{t) = E(0) 



k{2Af/ 2 t r< 2tl dm cos 3 ''I 2 (pm) 



7T/2M ^p(m,0) 



ttAAH 2 
7% - 1) 



(59) 



One can check that the decay of the thermal and macro- 
scopic kinetic energies, separately, is also quadratic in 
time. One also observes that, at the time of the density 
blowup, t — t c , the energy constitutes a finite and non- 
zero fraction of the initial energy, so the time t c is not 
special for the function E(t). 

As a simple example, consider the initial density 
p(m, 0) = po cos(pm) corresponding to Eqs. ([33)1 and 
(f36|) . Here the integration in Eq. ([59]) is elementary, and 
we obtain 



E(t) _ 7T 

p T l 7-1 



7T 

12 



47 



7T7 



2(7-1) 



(60) 

where poTqI is the characteristic energy scale. This E(t) 
dependence is shown in Fig. [TJ We now proceed to the 
analysis of the local structure of the flow near a singular- 
ity. 
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F. Local structure of the exact solutions near the 
singularity 

To analyze the local structure of the developing sin- 
gularity we consider the Taylor expansion of the inverse 
density it(m, t) = l/p(m, t) in a close vicinity of m — too 
at times close to t c . To calculate the m-derivatives it 
is convenient to write u(m,t) — u(to,0)[1 — tcf>(m)] 2 , 
where 4>(m) — /uy Ap(m, 0) cos(/zto) satisfies the condi- 
tions 0(to o ) — l/t c and 4>'(mo) = 0, see Eqs. ([25)1 and 
After some algebra we find 

u'(m ,t) a2 
u'(m ,0) 
u"(m ,t) 



u(m ,0) 
u"'(m ,t) 



-2t c 0"A + O(A 2 ), 

~2t c A [3ii'(m ,0)f 
+0(A 2 ), 



■u(m O) O)0" 



,.(4) 



(mo,*) 



u(m o ,0) 



= 6£ 



0(A), 



(61) 



where all the derivatives of <f> are evaluated at m = mo, 
and A = 1 — t/t c . The first non- vanishing m-dcrivative 
at t = t c is, therefore, of the fourth order, and we obtain 



u(m,t c ) 
u(m ,0) 



b 2 p 4 (m — m ) 4 , (i\m — m \ -C 1, (62) 



where we have introduced a positive dimensionless con- 
stant b — —t c (f>"(mo)/2fi 2 [recall that <j)(m) has a max- 
imum at to = too, so that </)"(m ) < 0]. Using the ex- 
pression for * c from Eq. (|2"9")l and the definition of one 
finds that b is independent of A and can be written as 



1 d 2 



where we have used [p(m, 0) cos(pm)]' (to = too) 
0. The latter relation implies (p'/p)(m = too) 
/j,tan(/iTO ), so we obtain 



1 + 2 tan TO 



p"(TO O ,0) 



fi 2 p(m ,0) 



Oil) 



(63) 



where the latter estimate assumes that the initial density 
p(m, 0) varies over a scale of order Equation (162)) 

shows that, at the singularity, u = 1/p vanishes faster 
than quadratically in to, as expected in general when a 
singularity is analytic. Going back to the Eulerian coor- 
dinates, x(m, t) — x(mo, t) — J mo u(m', t)dm' , we rewrite 
Eq. JB2J as 



u(x,t c ) 


5y/b(x — x c ) 


4/5 




u{x o ,0) ~ 


I 


5 


I 



1/5 



< 1, (64) 



where Xq = a; (too, 0), the spatial coordinate of the singu- 
larity x c is determined by Eq. (j2"5)) , and I — l/[pp(xo, 0)] 



is the inelastic cooling length scale. The validity condi- 
tion in x in Eq. (|64p corresponds to the validity condition 
in to in Eq. (|62| . Thus, at t — t c the density profile is 
singular in a vicinity of x = x c , and exhibits a power law 
with exponent 4/5: 



p(x o ,0) ~ \b\/b\x - x c \ 



4/5 



(65) 



This power-law singularity is integrable (that is, has a 
finite mass) and symmetric with respect to x c . For com- 
parison, the density singularity of a free flow, see e.g^, 
exhibits the exponent 2/3 instead of 4/5 . 

As the velocity itself is finite at singularity, the velocity 
gradient is of interest. We obtain 



dv 

dx 



p(m, t) 



x—x(m,t) 



dv 

dm 



2<t>(m) 
1 — t<j)(m) ' 



(66) 



where the last equality follows from the continuity equa- 
tion and definition of d>(m). As a result, 



2 



1 



cj)(m] 



— t ~ t c — t + bt c p 2 {m — mo) 2 , (67) 



up to higher order terms in t c — t and to — mo. Using 
Eqs. dnZJ) and we find that 



u(m,t c 



t c d x v(x,t c ) y 
This relation, combined with Eq 



dv(x, t c ) 
dx 



u(mo, 0) 
, (dH), yields 
I i 2/5 



5^(3 



(68) 



Note that while the exponent 2/5 of this power law is 
different from the exponent 4/5 of the power law for the 
density, the two power laws have the same region of va- 
lidity in x, see Eq. (|64|) . Though the velocity gradient 
diverges at the singularity, the velocity itself is continu- 
ous there and has a cusp ~ \x c — x\ 3 ^ 5 . This is in contrast 
with the "wave breaking" singularity of ordinary gas dy- 
namics, where the velocity becomes discontinuous at the 
point where the velocity gradient blows up. 

It is clear from the above that the local profiles of the 
density and velocity at t — t c do not depend on the details 
of the initial density p(m, 0). This is not so for the pres- 
sure for which two types of spatial behavior are possible. 
For the special case where p(m, 0) cos(/ito) is maximum 
at to = 0, and so too = 0, the density blows up at the ori- 
gin x — 0. As a result, the pressure versus x has a local 
maximum at x — 0, and the value of the maximum stays 
constant in time: p{x = 0, t) = 2A. The pressure profile 
in a vicinity of x = can be obtained from a Taylor ex- 
pansion of Eq. pp]): p(m, t)/A ~ 2 - p 2 m 2 + 0(p 4 m 4 ). 
At t = t c this leads to 



p(x,t c ) 1 / 5x 

p(0,t c ) ~ 2\lb 2 



2/5 



(?) 



2/5 



< 1. 



(69) 
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In the generic case, where the maximum of 
pirn, 0) cos(/xm) is not at zero, so that m ^ 0, 
the singularity develops at a point which is not special 
for the pressure. Here the Taylor expansion of Eq. (j2"0|) 
yields 

p{m)/A — 2 cos(/iTOo) — 2/i(m — mo) sin(/imo) 

+ O[p?(m-m f}. (70) 

In this case we find 

1/5 



corresponding to the condition p\m — mo| < 1 in Eq. 
([72")) . The asymptotes of F(z) are 



F(z) 



z, z 2 < 1, 

(hz/b 2 f'\ z 2 / 5 >l 



(77) 



The applicability condition (|76p is determined by the 
\z\ ^> 1 asymptote of F{z) and simplifies to 



P{x c ,t c ) 
P{x c ,t c ) 



1 — tan(/^m ) 



5(x - x c ) 



lb 2 



(71) 



x - x c (t) 



I 



1/5 



< 1 



that holds at (\x - x c \/l) 1/5 < 1. Note that the local 
x-dependence of the pressure at t = t c is very different 
from that of the density, or velocity gradient. First, the 
pressure remains finite at t = t c . Second, even though 
the pressure gradient diverges at x — x c , the divergence 
stems from a small correction term to a constant pres- 
sure, see Eq. (jB"9")) and (fTTj) . As we discuss later, this 
difference in behavior is crucial for understanding the 
physical nature of the singularity. 

Now let us investigate the local structure of the flow 
immediately before the singularity: at t c — t <C t c . The 
leading terms of the double Taylor expansion of the in- 
verse density u(m,t) = u(m, 0)[1 — t<p(m)] 2 in a vicinity 
of m — niQ and t = t c (that is, at fi\m — mo| <C 1 and 
A< 1) are the following: 



The same condition guarantees the validity of the power- 
law density profile (f65|) at t ~t c . 

Therefore, near the singularity the density has the fol- 
lowing self-similar form: 



1 



P(x,t) = 

p(x ,0) (l-t/tc) 2 



R 



x - x c (t) 



l{l-t/t c f/* 



(78) 



where R(z) = 1/[1 + bF 2 (z)] 2 . In the region correspond- 
ing toz<l one has R{z) ~ 1. That is, p(x,t) develops 
a plateau in a narrow region near x c (t) that we will call 
the inner region. The inner region shrinks as (t c — i) 5 / 2 
as t approaches t c : 



p{x,t) 



u(m, t) 
u(m o ,0) 



at 



\x - x c (t)\ 



1(1 - t/t c f/ 2 



< 1. (79) 



A 2 +2bA(j 2 (m-m a ) 2 +b 2 (i 4: {m-m ) 4: ; (72) At intermediate distances, or in the outer region, 



we recall that A = 1 — t/t c . Equation ([72]) can be written 
in a self-similar form: 



\x - X c (t)\ 



l(l~t/t c f/ 2 



2/5 



»1, 



X - X c {t) 



I 



1/5 



« 1 



it(m, t) 
u{m Q , 0) 



1 



jLt(m — mo) 



y/1 ~ t/tc 



(73) 



where U(y) = (1 + by 2 ) 2 . Now, x(m,t) 
J m u(m' ,t)dm' can be written as 



x c {t) 



the time-independent power law (p5]) builds up. At 
t = t c the power law rules in the whole region Qx — 

X c (t)|//) 1/5 « I- 

The development of the singular profile of the velocity 
gradient can be inferred by going back to Eq. (|67[) : 



x(m, t) — x c (t) 
I 



1 

tr 



5/2 



p(m — mo) 

\/i - Wc 



dv 
dx 



(t c -t)[l + bF 2 (z)Y 



[x - x c (t)} 
l{l-t/t c f/ 2 



(80) 



where ^(y) — Jq U(y')dy' . Evaluating this integral, we 
arrive at the self-similar form in the Eulcrian coordinates: 



u(x, t) 
u(x Ql 0) 



1 



1 + bF 2 



x c (t) 



l{l~t/t c f/ 2 



(74) 

where the function F(z) is denned implicitly by the fifth 
order polynomial equation 



Thus d x v develops a plateau d x v = —2/(t c — t) [cf. Eq. 
P2p ] in the inner region, while the power law described 
by Eq. (168j) sets in the outer region in the same way as 
the density power law. 

The development of the pressure profile as t approaches 
t c is different in the case of mo = and mo 7^ 0. In the 
former case we have 



b 2 F 5 (z) 2bF 3 (z) 



+ F(z) = 



(75) 



5 3 

that has a unique real solution. Equation (fT4"|) holds at 
x that satisfy the strong inequality 



p{x,t) 



1 



1 



Jills. F 2 



l(l-t/t c f/ 2 



(81) 



1/2 



F 



x c (t) 



1(1 - t/tc) 5 / 2 



< 1 



(76) 



that holds at (|x|/Z) 2 / 5 C 1. In the inner region, 

[\ x \/i(i-t/t c y>/ 2 ] 2 «i, 

p(x,t) 



1 



p(0,t) 



(1 - t/t c y 21 2 



(82) 
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while in the outer region the time-independent profile 
([55]) sets it. 

In the generic case of too ^ we obtain, for 
[\x-x c (t)\/l} 1/5 <$:!, 



p[x c (t),t] 



— tan(/xmo) yl — t/t c F 
In the inner region this yields 



.T — X c (t) 



p[a; c (*),*] 



1 — tan(/iTOo) 



a; c (t) 



?(l-iAc) 2 



(83) 



whereas the time-independent profile ([7T]) sets in the 
outer region. We observe that in the inner region, 
[\x — a: c (t)|//] 1 ^ 5 <C 1, the pressure is approximately con- 
stant. This suggests that the sound waves are the fastest 
physical process near the singularity. Before we consider 
the hierarchy of time scales in more detail, let us reiter- 
ate that the finite-time singularity described here is quite 
different from the free-flow singularity 24 where the den- 
sity blows up as (t c — t) , the plateau of the inner region 
shrinks with time as (t c — t) 3 / 2 , the power law tail in the 
outer region is ~ x 2 / 3 , and where the Lagrangian veloc- 
ity, rather than the Lagrangian acceleration, is constant 
in time. 



G. Time scale separation and isobaric scenario 

In general, there are three time scales that characterize 
the dynamics described by the nonlinear IGHD equations 
([5]) and ©: the sound travel time t sound ~ L/y/T, the 
cooling time t coo n ng ~ 1/ (ApVT) , and the inertial time 
tinertiai ~ L/v. Here p, T, and v are typical values of the 
fields, while L = L(t) is the characteristic spatial scale 
of the flow. The evolution of the hydrodynamic fields 
can produce time scale separation: a strong inequality 
between the time scales. Moreover, the hierarchy of the 
time scales can be different in different regions of space. 
Let us evaluate the time scales for our exact solutions 
(f2U|) - (|22p . Here there are only two independent time 
scales: the sound travel and the cooling time scales. This 
stems from the fact that the compressional heating and 
the inelastic cooling balance each other in the equation 
for the pressure, so that ti nert i a i ~ t cooling- Consider the 
inner region \x — x c \ < 1(1 — t/t c ) 5 ^ 2 . From the results 
of the previous subsection, yT ~ \/To(l — t/tc)i while 
L(t) ~ 1(1 - t/t c ) 5/2 . As ~ r ~ t c , we obtain 



tsound rs " 1 tc 



t 

l — 

tc 

-2 



3/2 



Now, using p ~ po(l — t/t c ) , we find 



cooling 



1 - — 

tc 



(84) 



(85) 



We observe that, except close to t c , t sound ~ t coo ( ing ~ t c . 
However, as the singularity is approached, the sound 
travel time in the inner region becomes much shorter 
than the cooling time. In this situation, the pressure 
in the inner region is expected to become constant^ as 
our solutions indeed show. 

To further elucidate this point, let us estimate the 
size of the spatial region at t = t c such that, within 
this region, the time scales obey the strong inequality 
tsound <C t cooling- Equations ([84]) and ([85]) show that, 
at a/1 — t/t c < 1, t soun d < t cooling in the inner region. 
At these times the size L(t) of the shrinking inner re- 
gion obeys the inequality (L/l) 1 ^ 5 -C 1. As we saw in the 
previous subsection, the shrinking inner region leaves be- 
hind stationary profiles of the fields. Therefore, at t = t c , 
the local time scales at some |x| = x Q <C I can be esti- 
mated by their values at those times when L(t) shrinks 
to the size xo. That is, at t = t c the time scale sep- 
aration tsound < tcooimg ~ tinertiai holds in the region 
(\x — x c \/l) 1 ^ <C 1 which is precisely the region where 
the pressure is approximately constant. Therefore, the 
density blowup, as featured by our exact solutions, lo- 
cally follows the isobaric scenario, previously suggested 
in the context of condensation processes developing in 
gases and plasmas that cool by their own radiation^ 3 ^. 
This fact has important consequences for the theory of 
clustering that will be explored elsewhere. 

We now proceed to the derivation of an additional so- 
lution that would allow us to demonstrate two important 
features of the developing density singularity: its locality 
and its possible coexistence with shock singularities of 
ordinary gas dynamics. 



IV. A PISTON MOVING INTO A GRANULAR 
GAS AT REST: A BLOWUP WITH A SHOCK 

The family of exact solutions, describing the finite time 
density blowup and reported in Sect ion InTl have a special 
value of the total gas mass, n/fi = (V2tt~/) / A. Here we 
relax this requirement by constructing an exact solution 
that can have an arbitrarily large mass. In this solution 
both the finite-time density blowup and an "ordinary" 
shock discontinuity are present. 

First, we note that formation of a density singularity 
in hydrodynamics is a local process. The set of Eqs. ([TO]) - 
(|llj) is hyperbolic and has a finite speed of propagation 
of information. Therefore, a finite-time density blowup, 
developing at a point with a finite x, cannot be affected 
by a change in the initial conditions sufficiently far away 
(provided the velocity is finite there). In particular, this 
is true for initial density variations that change the total 
mass of the gas, possibly making it infinite. The solution 
that we are going to present here illustrates this point, 
as it has an arbitrary density distribution sufficiently far 
from the developing singularity. 

The solution also illustrates another point that is ab- 
sent in the solutions reported in Section [HI! appearance 
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of shocks. At A = 0, Eqs. (p~0|) - (p~T|) become the equations 
of classical ideal gas dynamics which produce shocks: ini- 
tially smooth hydrodynamic fields develop shock discon- 
tinuities in a finite tim o 26 ' 27 . Shocks can also appear at 
A > 0, the case of our interest. The following argument 
can be helpful in elucidating the comparative role of the 
two types of singularities: the density blowup and the 
shock. Let the initial conditions be fixed. Then, as A 
goes down, the development of a density blow up will 
be delayed in time (the delay time becoming infinite as 
A — > 0). On the other hand, the time of shock forma- 
tion must obviously approach a finite limit as A — > 0. 
Therefore, for sufficiently small A the shock formation 
will typically precede the density blowup. In the exact 
solution that we are going to present, the shock is present 
in the solution from the very beginning. 

We choose for this solution (an extended version of) 
a basic setting of ideal one-dimensional gas dynamics: 
at t = a piston starts moving at a constant speed Vq 
into an undriven granular gas at rest. Because such a 
gas has a zero temperature everywhere, the initial state 
of the gas is uniquely characterized by the initial den- 
sity profile, say po(x). It is convenient to go over to 
the frame moving with the piston, where the piston rests 
at x = 0. There one needs to solve Eqs. §5§ and (J6]) 
with the initial conditions p(x, 0) = po{%), v(x, 0) = — Vq, 
T(x, 0) = and the boundary conditions v(x = 0, t) = 
and v(x — +oo,t) = — vq. At t = the gas hits the 
piston wall at x = 0, and a shock forms instantaneously 
and starts propagating into the gas. Each of the hydrody- 
namic fields experiences a discontinuity at the shock front 
Xo(t) that obeys Xq(0) = 0. The solution at x > Xo(t) 
is of course p(x,t) = po(x + Wot), v(x,t) = —vq and 
T(x,t) = 0, while at < x < Xo(t) non-trivial distribu- 
tions of the hydrodynamic fields develop. We show below 
that a special choice of po{x) yields a solution that, at 
x < xq (t) , is of the same type as that described in Section 

nm 

The jump conditions at the shock front are provided by 
the same Rankine-Hugoniot conditions as in the ordinary 
gas 2 ^. Indeed, these conditions can be obtained by con- 
sidering the equations of mass, momentum and energy 
balance in an infinitesimal volume (xo(t) — e, Xo(t) + e), 
in the limit of e — ► 0. While the mass and momentum 
balances are exactly the same as in the ordinary gas, the 
inelastic loss correction to the energy balance is propor- 
tional to f*°{$** p 2 T 3 l 2 dx (see subsection MM > an d it 
vanishes at e — ► despite the discontinuity of the in- 
tegrand. The Rankine-Hugoniot conditions state that, 
in the frame moving with the piston, the mass, momen- 
tum and energy fluxes should be continuous through the 
shock 2 ^: 



p 2 (v 2 - Xq) = 
p 2 (v 2 - X ) 2 + 

p 2 (v 2 - X ) 3 + 



-pi(v + x ), (86) 
p 2 T 2 = pi(u + x ) 2 , (87) 
27P2(«2 - x )T 2 
7-1 



~Pi(vq + Xq) 3 



where the subscripts 1 and 2 stand for the upstream and 
downstream values of the fields, respectively. These con- 
ditions yield 



dx = (7 - j>o (7 + 1>2 
dt ~ 2 2 



(89) 



(7 + rr (7-1) fa) + ^2) , on x 
Pi= 7 _! » T 2 = 2 • (90) 

Now let us go over to the Lagrangian coordinates, see 
Eq. (fl4)) . We denote the Lagrangian coordinate of the 
shock front by mf(t): 

/.K (t)-<5 

m/(t) = lim / p(x',t)dx', (91) 

where 5 > 0. Differentiating Eq. (f9"Tj) over t we find 

' dxo 



(92) 



drrif ( dxo 



where the second equality holds by virtue of Eq. 
Using this result together with Eqs. ([89]) and ([90]) . we 
obtain 



drrif 
dt 



'(7 + l)P2 



2«i 



u 2 



7' 



1 



7 + 1 



til, 



V2 



-V 



2p 2 Ui 

7 + 1 



(93) 



where the inverse density tti : 2 = l/pi.2 is introduced and 
T is expressed via p and u. Now we show that it is possi- 
ble to choose the initial density profile po{m) [or, equiv- 
alently, po{x)] so that the solution of Eqs. (HHJ) and fTT]) 
at < to < mf (t) is a constant acceleration solution pre- 
sented in Sect ion ITTT1 In the upstream region, m > m/(t), 
the gas is undisturbed by the shock, and the solution is 
u(m,t) = uo(m), v(m,t) = —vq and p(m,t) = 0. In the 
downstream region < 777 < rrif(t) the hydrodynamic 
fields at t > are 



p(m) = 2 A cos(pm), 



u(m,t) = /(m) — pt^/ Acos(pm) , 
v(m,t) = — 2/i / /(to')y A cos(pm')dm' 



2 Apt sinfam), 



(94) 
(95) 

(96) 



where /(to) is a yet unknown function. By construction, 
the upstream and downstream solutions obey the govern- 
ing equations, and what is left is to impose the boundary 
conditions (|9"3"|) at the shock front. The first condition 
becomes 



drrif A(j + 1) cos[prrif(t)] 



dt 



u [m f (t)] 



(97) 
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Once uo(m) is known, this equation, together with the 
initial condition mf(0) = 0, determines the shock coor- 
dinate versus time, rrif = m-f(t), via its inverse function: 
t = to(mf), where 



to(m) 



ito(m') 



o y A(7 + 1) cos(/xm') 



-dm'. (98) 



Using the inverse function to(rn), we demand the second 
and third conditions in Eq. (|93[) . 



/(m) - /i*o(m)VAco S ( M m) = J (7 l |" 1 ° (m) 

V 7+1 



(99) 



-2/i / f(m')y/Acos(pm')dm' + 2^4/zio (to) sin(/im) 
Jo 



' 4Auo (?n) cos(pm) 



7+1 



(100) 



on the interval < m < 7r/(2/i). Using Eqs. (f9"8")) and 
(f9"9"| . we can express /(m) via rto(m): 



/(m) = /i / 
Jo 



uo(m') cos(/Um) , 
(7 + I) cos(fim') 



(7 - l)wo(m) 



7+1 



(101) 



Note that this relation does not include the parameter A. 
We now substitute Eq. (|101[) into Eq. (| 1 00|) and arrive at 
a closed equation for uo(m): 



-2p 2 I dm cos(/im') / cim" 



uo(to") 



(7 + 1) cos(/im") 



, / (7 — l)ito(m') cos(//m') 



-2u / dm , , 

7o V 7+1 



+2/i sin(/xm) 



uo(m') 



o y (7 + 1) cos(pm') 



dm' 



vo / 4wo(m) cos(/xm) 

3 + y 7 + 1 ' 



(102) 



This cumbersome equation is a linear integral equation 
for the function \J i*o(m), and it is soluble analytically. 
Changing the order of integration in the first term of the 
equation, we rewrite the first term as 



-2fi 2 / dm"J-. -Hi—/ — - 

io V 7+1 cos /zm" 



coa(fj,m')dm' 



-2usin(//m) / dm" \ . 

io V (7 + l)cos(/zm") 



+2/x / dm" s'm(pm")^ 



uo(m") 
- 1) cos( 

wo (m") 



(7 + 1) cos(/xm") 



This brings a partial cancelation of terms in Eq. (|102[) . 
and we obtain a simpler equation 



2u / dm' sin( um')\ -, -H- — — 

Jo V (7 + 1) cos(Aim') 



-2/X4 



' (7-1) 
7+1 



dm! \J uo(m') cos(pm') 



/ i«o(to) cos(pm) 
M V 7+1 ' 

Now let us introduce an auxiliary function 



g(m) 



I uo(m) cos(pm) 



7+1 



that obeys, on the interval < m < 7r/(2/i), a linear 
integral equation: 



g(m') tan(/^m/) — \J ^ — 1 



dm! 



vo 



The solution for g(m) is elementary: 

w exp(-V7~ 
g(m) = ■= , 

y/A cos(pm) 

so the result for uo(m) is 

v o(l + 1) ex P ( — 2 V7 - 1m to ) 



u (m) 



4A cos 3 (/Ltm) 



+ 5( m ) 
(103) 

(104) 
(105) 



To complete the formal construction of the solution, we 
present the initial gas density in the Eulerian coordinates, 
po(x), in a parametric form: 



Po(m) — po cos 3 (jim) exp ^2^/7 — 1 prn^j , 

dm' 
Po(m') ' 



(106) 



where po — 4A/[(7 + 1)wq]. The graph of po(a^) is shown 
in Fig. [5l 

Now we use Eq. (|105l) to calculate the inverse func- 
tion t — to(mf) from Eq. (f98|) that determines the shock 
motion law in the Lagrangian coordinates, mf = mf(t): 



t (m) = t*$(p,m), $(z) = 



exp (— V7 — 1 z ) dz' 



(107) 

where we have introduced the characteristic inelastic 
cooling time scale i» = 2Z/[(7 + l)i>o], and I = l/(pp ) 
is the inelastic cooling length scale. As $>(z) diverges at 
z = tt/2, mf(t) satisfies, at any finite time, the double 
inequality < m/(t) < n/(2p). We reiterate that, in the 
upstream region m > m/(t), the gas is unperturbed by 
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vo(t/t*) s'm(pm) 



(108) 



FIG. 5: The initial density po(x) (in units of po) versus the 
Eulerian coordinate (in units of I), see Eq. (|106[) , in the 
problem of a piston moving into a granular gas at rest for 
7 = 2. 




FIG. 6: (Color online) The exact solution (| 108|) . in the La- 
grangian coordinates, of the problem of a piston moving into 
a granular gas at rest. Shown are the inverse density (a), ve- 
locity (b), pressure (c) and temperature (d) at rescaled times 
(from left to right) tvo/l = 0.2, 0.4, 0.6 and 0.65 as a func- 
tion of m measured in units of 1/fj, (the solid lines). The 
dashed lines are the initial conditions: the initial density from 
Eq. (|106p . the constant (rescaled) velocity —1, and zero tem- 
perature and pressure. The density, velocity, pressure and 
temperature are rescaled to po, vo, poVo and Do, respectively, 
and 7 = 2. The density blows up at the piston, m — 0, at the 
rescaled time t c vo/l = t*Vo/l = 2/3. The rescaled Lagrangian 
coordinate of the shock at the blowup time is 1.08031 . . .. 



the shock. The final form of the solution in the down- 
stream region, < to < m/(t), is 



p(m,t) = 
p(m,t) = 

v{m,t) = 



(7 + 1 )PoVq cos(/uto) 



cos(pm) [&(fim) + V7 - l$'(£tm) - t/U\ 

—vq / dz cos z 
Jo 



$(z) + v/ 7 - l$'(z) 



This solution is shown in Fig. [5] Note that, at fixed po 
and vq, the characteristic spatial and temporal scales of 
the solution behave as 1/A. For this solution, the gas 
density blows up in a finite time at the point m where 
the function $(/im) + y/j — l<$>'{p,m) has its minimum. 
One can easily see that, for 7 < 2, this function is mono- 
tone increasing with to. As a consequence, the density 
blows up at the piston (that is, at m = 0), and this 
happens at t = t c = ^7 —It*. At this time, the shock 
front location to* = mf(t») is described by the relation 
$(/jm,) = y/j — 1 which yields fim, = 1.08031 . . . for 
7 = 2 (in 2d) and 0.87915 ... for 7 = 5/3 (in 3d). That 
is, by the time t = t c when the density blows up at 
the piston, the shock has traveled only a finite distance 
(both in the to-, and in the x-space) from the piston. It is 
obvious, therefore, that our solution allows an arbitrary 
modification of the density profile po(x) at sufficiently 
large x that are unreachable for the shock. This clearly 
shows that initial states with an arbitrary large mass of 
gas can develop a finite-time density blowup. 



V. NUMERICAL SOLUTIONS 

We confirmed the exact solutions, presented in Figs. 1- 
4 and 6, by solving numerically the IGHD equations (fT6j) 
and (fT7|) in the Lagrangian coordinates. In each case 
we used the initial and boundary conditions, provided 
by the exact solutions themselves. We employed the 
classical artificial viscosity, staggered grid scheme of von 
Neumann and Richtmyer 31 . The number of cells (grid 
points) used varied between 1000 and 2000 and showed 
numerical convergence until close to the singularity. As 
the solution developed a density blowup, the simulation 
had to be terminated at a very high but finite maximum 
density: usually at about 10 7 poj where po is the initial 
density. We also enforced the simulations to stop when 
the density jump between two adjacent cells exceeded a 
prescribed value, usually 50 percent. (When such jumps 
develop near the singularity, the accuracy of the simu- 
lation degrades and can be restored only by a rezoning 
algorithm that was not employed.) One example of the 
numerical solution is shown in Fig. [T] alongside with the 
analytical solution, and very good agreement is observed. 
Very good agreement was also obtained for the solutions 
shown in Figs. 3, 4 and 6 (not shown). The numerical 
solutions inevitably add some effective noise to the sys- 
tem because of the spatial and temporal discretization 
and round-off errors. The fact that the analytical so- 
lutions are accurately reproduced numerically confirms 
their stability with respect to small one-dimensional per- 
turbations. 

We also performed extensive numerical simulations 
with Eqs. (fT5|) and ((TTJ) for a variety of initial and bound- 
ary conditions that did not correspond to any of the 
exact solutions. As already briefly reported in Ref. 22 , 
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FIG. 7: (Color online) An example of numerical solution of 
Eqs. (|16|) and 1)170 . The initial conditions are p(m,0) = po, 
T(m,0) = T and v(m,0) = -v tanh [m/(p X)]. (a) The 
inverse density profiles at rescaled times 0, 2.06, 3.76, 5.36 
and 6.26. (b) The density p(x,t) in the logarithmic scale at 
times t = 2, 4, 5, 5.5 and 7 (from the bottom to the top, 
respectively) as a function of x. (c) The time history of the 
density (the solid line) and pressure (the dashed line), both 
in the logarithmic scale, at the blowup point x — 0. The inset 
depicts p _1//2 (s = 0,t) in the logarithmic scale versus time in 
a vicinity of t = t c . (d) The rescaled density p(x,t)(t c — t) 2 
as a function of the rescaled coordinate x(t c — 1)~ 5 ^ 2 at times 
t = 6.7, 6.9, 7, 7.05 and 7.07 (from the top to the bottom, 
respectively) in a vicinity of the singularity. The profiles at 
different times almost coincide demonstrating the local self- 
similarity of the blowup. (e) The log-log plot of the density 
p(x,t) versus x (solid lines) for the same times as in (c) (from 
the bottom to the top, respectively). The dashed line shows 
a power law p ~ x~ 4 ^ 5 , to guide the eye. The Lagrangian and 
Eulerian coordinates, time and the hydrodynamic fields are 
made dimensionless, as explained in Section [V] The parame- 
ters are A = poAA = 0.5, Vq — \/To and 7 = 2. 



these simulations show that, for generic initial condi- 
tions, a finite-time density blowup always occurs. Re- 
markably, the numerical solutions exhibit, close to the 
singularity, the same local scaling behaviors in space 
and in time as those exhibited by our exact solutions 
and presented in subsection IIII F) One series of sim- 
ulations dealt with a symmetric inflow of an initially 
uniform gas, p(x,0) — po and T(x,0) = Tq, from "in- 
finity": v(x,0) — — vq tanh(x/A). Here it is convenient 
to rescale the x-coordinate by A, time by \/\/Tq, the 
density by po, the velocity by y/To, the temperature 
by T , the pressure by poT and the Lagrangian mass 
coordinate m by p X. After this rescaling the govern- 
ing equations Eqs. (fToj) and (fTT)) keep their form, ex- 
cept that A becomes rescaled: A = pqAX. The rescaled 



initial conditions become p(x,0) = p(x,0) = 1 and 
v(x,0) — — (vo/-\/To) tanh.(a;). The numerical solutions 
were obtained on the (rescaled) interval |x| < 10 that 
corresponds, at t = 0, to the (rescaled) Lagrangian in- 
terval \m\ < 10. The boundary conditions are v(x — 
±10, t) = v(x = ±10,0) = T(V\/Jo)tanhl0. Because 
of the symmetry of the problem with respect to x — we 
actually solved it on the half-interval < x < 10, replac- 
ing the boundary condition at x = —10 by the condition 
v(0,t) = 0. 

Here we present one typical example of such a simu- 
lation, briefly mentioned in Ref^. The parameters are 
A = 0.5, v = V% and 7 = 2. The (one half of the) 
simulated flow is shown in Fig. [7] Panel (a) provides a 
general view of (one half of) the system. The gas inflow 
creates a compression in the central region. A compres- 
sion wave propagates outwards and steepens. This steep- 
ening would lead to a wave breaking singularity, but the 
numerical scheme resolves this singularity, by means of 
the artificial viscosity, as a shock wave that continues 
propagating outward. At the same time the gas density 
at the origin continues growing and ultimately blows up. 
The density growth in a vicinity of the origin is presented 
in panel (b). Panel (c) focuses on the density and pres- 
sure history at the origin. While the density blows up at 
t = t c , the pressure hardly varies there so the isobaric sce- 
nario holds. The inset of panel (c) verifies that, close to 
t c , the density blowup proceeds as (t c — t)~ 2 . Panels (d) 
and (e) add more tests of the spatial and temporal scal- 
ing behavior near the singularity. Panel (d) shows a plot 
of the rescaled density p{x, t)(t c — t) 2 versus the rescaled 
coordinate x(t c — i)~ 5 / 2 . The collapse of the curves at 
different times proves the local self-similarity (note that 
the gas density at the origin varies, for these times, by 
four orders of magnitude) . Finally, panel (e) verifies the 
presence of the inner and outer regions, described by our 
theory, see subsection IIII Fl The density plateau, whose 
size shrinks as ~ (t c — t) 5 ^ 2 , represents the inner region, 
while in the outer region a time-independent density pro- 



as pre- 



file forms with a power-law behavior p ~ x 4 ^ 5 
dieted by our exact solutions. 

The same universal features of the singularity were also 
observed when starting from small-amplitude sinusoidal 
density or velocity perturbations around a homogeneous 
state. 



VI. NON-IDEAL EFFECTS NEAR THE 
SINGULARITY 

Having found the exact nonlinear solutions of the 
IGHD equations and (O , we are in the position to test 
the assumptions behind these equations and establish the 
domain of validity of the solutions (as accurate leading- 
order descriptions) in the presence of "non-ideal" effects. 
First, the validity of the Navier-Stokes hydrodynamics, 
Eqs. ©-(HI), demands that the Knudsen number lf ree /L 
remain much smaller than unity. The smallest hydrody- 
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namic length scale of the solution can be estimated as 
L(t) ~ l(l-t/t c ) 5 / 2 , see subsections UlTF] and HHni The 
mean free path lf ree ~ l/(pa d ~ 1 ) ~ (1 — t/t c ) 2 / (p a d ~ 1 ). 
As I ~ l/(Ap ) and A ~ (1 - r)CT d_1 , we obtain 



^ /re 



l-r 



VI - */*c 



(109) 



At i < t c the Knudsen number is small by the assump- 
tion 1 — r <C 1. However, as t approaches i c) the Knudsen 
number grows indefinitely which invalidates the hydrody- 
namics. We shall see, however, that one of the assump- 
tions of the IGHD breaks down even earlier. 

Now consider the ratio of the viscous stress term to 
the pressure gradient term in the momentum equation 
([3]). This ratio can be estimated as (voVTd x v) / (pT) = 
Vodmv/vT. We first estimate all the ratios in the inner 
region, see subsection MI Fl As both d m v and yT van- 
ish linearly as t — > t c , their ratio depends on time only 
weakly, and can be estimated by its value at t = 0. As 
\d m v(m,0)\ ~ Ay/T (see subsection IIII Al) . we find 



(110) 



i>od x (VTd x v) 



d x { P T) 



that is the viscous stress is always negligible in our solu- 
tions as long as 1 — r <C 1. 

The same estimate (up to the sign) holds for the ratio 
of the viscous heating term to the compressional heating 
term in Eq. (gj): 



Mi - l)VT(d x v) 2 / P 



( 7 - l)Td x v 



1 - r < 1 



(111) 



In contrast to the viscous terms, the heat conduction 
term in Eq. ((4|), which is initially small in our solutions, 
does become important near the singularity. We estimate 
the ratio of the heat conduction term to the inelastic 
cooling term as follows: 



K d x (VTd x T)/p 



ApT 3 / 2 



Kq 



l-r 



Ap 2 L 2 1 - t/t c 



(112) 



The same estimate is obtained for the ratio of the heat 
conduction term to the compressional heating term. This 
ratio becomes of order unity at 1 — i/t c ~ 1 — rCl. At 
this time lf ree /L ~ y/l — r -C 1, see Eq. (|109[) . so the 
IGHD equations break down while the full Navier-Stokes 
hydrodynamics, Eqs. (J2|)-(j4|), is still valid. 

In its turn, the dilute gas assumption, pa d <C 1, breaks 
down, and excluded particle volume effects become im- 
portant, at 1 — t/t c ~ y/ poa d <C 1. The relative role 
of the heat conduction and excluded particle volume ef- 
fects in the breakdown of our analytic solutions near 
the attempted singularity is determined by the compe- 
tition between the small parameters y p$o d and 1 — r. 
When y/p (T d > 1 — r, our solutions remain valid until 



the density becomes comparable with the (fraction of) 
close packing density. This happens at 1 — t/t c ~ \j pa d , 
so L va nd ~ l(po<7 d ) 5 / 4 . At moderate densities one can 
use, in numerical solutions, a half-empiric equation of 
state due to Carnahan and Starling 32 , and half-empiric 
transport coefficients obtained for granular gases^ 3 - in the 
spirit of Enskog kinetic theory 3 ^. 

When y Po<^ d "C 1 — r the heat conduction becomes 
important when the gas is still dilute. This happens at 
time 1 — t/t c ~ 1 — r so that L va ud ~ 1(1 — r) 5 / 2 . As long 
as y/l — r <C 1 , the Knudsen number is still small at that 
time, and the complete Navier-Stokes hydrodynamics is 
still applicable. 

Therefore, as the blowup time t c is approached, either 
the heat conduction, or the excluded particle volume ef- 
fects become important and invalidate our theory in the 
inner region. In this sense, our solutions describe an in- 
termediate asymptotic regime of formation of a dense 
cluster of particles. Importantly, in the outer region our 
solutions remain valid until much later times. Indeed, as 
the inner region shrinks, it leaves behind stationary pro- 
files of the fields, see subsection IIII Fl The ratios of the 
different terms governing the IGHD validity are given, 
at some x from the outer region, by the corresponding 
ratios in the inner region estimated at the time when 
L(t) ~ \x — x c \. As a result, at some time close to t c 
(see above), our solutions break down in the inner re- 
gion which size at that time is L va nd, while in the outer 
region, \x — x c \ 3> L va ud, our solutions continue to hold. 

The fact that the IGHD description breaks down only 
locally, in a small region of space, can be exploited for 
derivation of an effective description of the clustering dy- 
namics beyond the singularity time. In this description 
close-packed clusters appear as finite-mass point-like sin- 
gularities of the density^ 2 .. This effective description is 
similar in spirit to the description of shocks (that ac- 
tually have finite widths) as discontinuities in ordinary 
ideal gas dynamics. 



VII. SUMMARY 

Let us briefly summarize the main results of this work. 
We introduced "ideal granular hydrodynamics" (IGHD) 
equations that provide a simple but informative descrip- 
tion of non-stationary large-scale flows of dilute granular 
gases with nearly elastic collisions between the particles. 
We employed the IGHD equations to investigate analyt- 
ically and numerically the paradigmatic phenomenon of 
particle clustering in freely cooling granular gases. We 
believe that the IGHD will provide a useful framework 
for a host of other problems involving large-scale flows of 
driven granular gases. 

We focused on a one-dimensional granular hydrody- 
namic flow, characteristic of an idealized channel geom- 
etry, and derived a family of exact nonlinear and non- 
stationary analytic solutions of the IGHD that describe a 
finite-time blowup of the gas density. The derivation was 
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made possible by a transformation of the hydrodynamic 
equations to Lagrangian coordinates. The exact solutions 
are characterized by a constant in time accelerations in 
the Lagrangian frame, and they are not self-similar. We 
investigated the local structure of the flow near the singu- 
larity and determined local spatial and temporal scaling 
laws for the hydrodynamic quantities in question. We 
also found an instructive soluble case of the problem of 
a piston entering, at a constant speed, a granular gas 
at rest. Here a density singularity, developing at the 
piston, coexists with a shock wave that is located else- 
where. Besides demonstrating the presence of the two 
different types of singularities in the same system, this 
solution allows an arbitrary density profile at large dis- 
tances, showing that the density blowup is a local pro- 
cess. In all of these solutions the developing singularity 
obeys locally the isobaric scenario^^: the gas pressure 
remains approximately uniform in space and constant in 
time in a close vicinity of the blowup point. This finding 
has important consequences for the theory of clustering. 
Indeed, by imposing, in the zeroth order of theory, the 
isobaricity condition on the hydrodynamic equations, one 
arrives at a powerful reduced description of the nonlinear 
clustering process. This description is valid for a much 
broader class of initial conditions than those giving rise 
to the exact solutions reported in the present work. The 
corresponding results will be presented elsewhere. 

Our numerical solutions of the IGHD equations accu- 
rately reproduce the analytic solutions, thus confirming 
their stability with respect to small one-dimensional per- 
turbations. Furthermore, numerical simulations with a 
variety of initial and boundary conditions showed that 
that the local scaling laws near the singularity are uni- 
versal, that is independent of details of the initial and 
boundary conditions. 

We also analyzed additional physical effects, neglected 
in IGHD, which become important in a narrow spatial 
region near the attempted singularity. Depending on the 
parameters, either excluded particle volume effects, or 
heat conduction invalidate our solutions there. In the 
former case, the density growth should directly cross over 



to the formation of close-packed clusters of particles. In 
the latter case, the final outcome of the nonlinear density 
growth is yet unknown. The future work should find out 
whether the heat conduction arrests the density blowup 
or only modifies the singularity law. In any case, our 
exact solutions can be viewed as instructive intermediate 
asymptotics, describing a broad class of strongly nonlin- 
ear clustering flows of granular gases. No less important, 
as breakdown of these solutions occurs only in the narrow 
regions around the density peaks, it is possible to con- 
tinue the ideal solutions beyond the singularity by intro- 
ducing into the theory finite-mass point-like singularities 
of the density. This procedure yields a simple effective 
description of granular gases with embedded close-packed 
clusters, in much the same way as the ideal gas dynamics 
describes a gas flow with shock discontinuities^. 

Very recently, one of our exact solutions has been suc- 
cessfully tested against MD simulations of the dynamics 
of a very dilute gas of nearly elastically colliding identical 
hard disks in a very long and narrow 2d channel. The 
results will be presented elsewhere. 

Finally, we note that the power law of the density 
blowup near the singularity, ~ (t c — <)~ 2 , is the same as 
the one recently found in a different setting: that of an 
initially thermalized granular gas, freely cooling and col- 
lapsing under gravity^. In each of the two problems the 
momentum equation is "fast" , while the energy equation 
is "slow" in a small vicinity of the collapse. This suggests 
common physics behind the two collapse phenomena and 
demands further investigation. 
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